Comparative transcriptome analysis of skin color-associated genes in leopard coral grouper (Plectropomus leopardus)

Background The leopard coral grouper (Plectropomus leopardus) is an important economic species in East Asia-Pacific countries. To meet the market demand, leopard coral grouper is facing overfishing and their population is rapidly declining. With the improvement of the artificial propagation technique, the leopard coral grouper has been successfully cultured by Fisheries Research Institute in Taiwan. However, the skin color of farmed individuals is often lacking bright redness. As such, the market price of farmed individuals is lower than wild-type. Results To understand the genetic mechanisms of skin coloration in leopard coral grouper, we compared leopard coral grouper with different skin colors through transcriptome analysis. Six cDNA libraries generated from wild-caught leopard coral grouper with different skin colors were characterized by using the Illumina platform. Reference-guided de novo transcriptome data of leopard coral grouper obtained 24,700 transcripts, and 1,089 differentially expressed genes (DEGs) were found between red and brown skin color individuals. The results showed that nine candidate DEGs (epha2, sema6d, acsl4, slc7a5, hipk1, nol6, timp2, slc25a42, and kdf1) significantly associated with skin color were detected by using comparative transcriptome analysis and quantitative real-time polymerase chain reaction (qRT-PCR). Conclusions The findings may provide genetic information for further skin color research, and to boost the market price of farmed leopard coral grouper by selective breeding. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-09091-6.

showed that overfishing and the destruction of spawning aggregations significantly declined leopard coral grouper population in the Philippines and Australia [6][7][8]. Even though the Red Book of the International Union for Conservation of Nature (IUCN) currently lists leopard coral grouper as the least concern [9], it is still essential to establish fishery management and conservation strategies to prevent leopard coral grouper from becoming a species of high-risk level in the future.
Coloration plays an important role in many creatures, which is associated with thermoregulation, social communication, predator avoidance, camouflage, protection from radiation, and selective mating [10,11]. Previous studies indicated that coloration of teleost fish is determined by six types of pigment cells, including melanophores (black and dark brown), xanthophores (yellow), erythrophores (red and orange), iridophores (reflective), leucophores (white), and cyanophores (blue) [12][13][14]. Studies also suggested that cAMP (cyclic adenosine monophosphate), MAPK (mitogen-activated protein kinase), PI3k/Akt (phosphatidylinositol 3-kinase-Akt), and Wnt (wingless-type MMTV integration site) signaling pathways are involved in melanin-synthesis related pathways in fish [15][16][17][18]. The leopard coral grouper has bright red skin color, which is different from the other species in the same genus. Due to bright red being associated with good luck and happiness in Chinese culture [19], the net price of bright red leopard coral groupers is approximately US$7 per kilogram higher than the brown individuals in the Asian market [5]. As such, for leopard coral grouper, skin color is an important economic trait. Although artificial breeding has been successfully established for leopard coral grouper by Fisheries Research Institute in Taiwan, mostly the skin color of farmed individuals is gray and brown. Several factors, such as habitats, physiological response, food, and genetics, were associated with skin coloration [20][21][22]. Feed additives and breeding approach have been proposed in order to improve their skin color performance [23,24]. Previous studies showed that the feed additives of carotenoids contain astaxanthin [4,25,26], enhancing the skin color of fish appearing red pigmentation. Besides, gene regulation associated with skin color is also a critical factor. Candidate carotenoid-related genes, including BCO2, LRP11, ANGPTLs, ALAS, PDIL, MED12, SOX, FAX, FATP, SCD, and LDLRA, were involved in carotenoid metabolism and significantly linked to the skin coloration of leopard coral grouper [5,22,27,28]. Several miRNAs including miR-215, miR-2188, miR-194, miR-122, and novel-m0118 may also play a potential role in regulating skin color in the red-colored leopard coral groupers [29].
Transcriptome analysis has been widely used to reveal gene expression affecting the skin coloration in nonmodel organisms, such as common carp (Cyprinus carpio) [17,40,41], crucian carp (Carassius auratus) [42], crimson snapper (Lutjanus erythropterus) [43], and catfish-like loach (Triplophysa siluroides) [44]. However, the molecular mechanism of skin coloration is still poorly understood in the Serranidae family. In this study, we characterized the transcriptome sequences of leopard coral grouper with different skin colors (red and brown) collected from Penghu Sea in Taiwan (Fig. 1). We aimed to (i) conduct transcriptome analysis of leopard coral groupers with different skin colors; (ii) validate differentially expressed genes (DEGs) associated with different skin colors of leopard coral groupers from other populations. The purpose of this study was to reveal novel information on candidate DEGs related to skin coloration between red and brown leopard coral groupers. Ultimately, the findings could provide useful genetic information for the artificial breeding of leopard coral grouper with different skin colors.

Illumina sequencing and genome-guided de novo transcriptome assembly
The transcriptome of six leopard coral groupers, including three red individuals and three brown individuals, were sequenced, respectively ( Clean reads with Q20 and Q30 were more than 92%, and the average GC content was 47.5%. The results showed that high-quality reads were used for transcriptome assembly in this study ( Table 1).
The final contig assembly produced 24,700 transcripts with an N50 length of 2,496 bp, an average contig of 1,776 bp, a medium length of 1,230 bp, and a GC content of 46.7% using CD-HIT-EST. The length distribution of all the transcripts was shown in Fig. 2. The longest and shortest contigs were 65,925 bp and 306 bp, respectively. There were 14,860 transcripts longer than 1,000 bp, accounting for 60% of complete data (Fig. 2).

Functional annotation
To annotate these assembled transcripts, which were based on UniProt Knowledgebase (UniProtKB), Cluster of Orthologous Groups database (COG), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genome Orthology (KO), including 14,171 genes, 5,492 genes, 6,411 genes, 8,480 genes accounting for a total of 24,700 genes, were obtained, respectively. 3,049 common genes were annotated by UniProtKB, COG, GO, and KO (Fig. 3). In the COG, the predominant categories were general function prediction only (20%), followed by cell cycle control, cell division, chromosome partitioning (17%), and function unknown (12%) (Additional file 1: Fig. S1). GO was divided into three major function categories: molecular functions (52%), cellular components (8%), and biological processes (40%) (Additional file 1: Fig. S2). KO classification obtained seven categories including 47 pathways. Lipid metabolism (108 genes) and transport and catabolism (128 genes) were the most abundant term in the categories of metabolism and cellular processes (Additional file 1: Fig. S3).

Genetic relationship between different skin color fishes
The principal component analysis (PCA) plot showed that the red (R1, R2, and R3) and brown individuals (B1, B2, and B3) were clearly separated with main principal component (PC) scores as follows: PC1 = 34.01% and PC2 = 19.55% (Fig. 4). Heatmaps based on the trimmed mean of M-values (TMM) expression of top 230 DEGs (115 up regulation and 115 down regulation) were generated from six transcriptome sequencing samples. Results showed that the transcripts count of red individuals (R1, R2, and R3) was different from that of the brown individuals (B1, B2, and B3) in Fig. 5 and Additional file 2: Table S1.

Differential expression genes
To analyze the DEGs between the red group (R1, R2, and R3) and the brown group (B1, B2, and B3). There were 1,089 DEGs (613 up-regulated and 476 downregulated DEGs) obtained based on an absolute value of log 2 fold change (FC) ≥ 2 and false discovery rate (FDR) < 0.05 (Additional file 1: Fig. S4). For the GO analysis, up-regulated and down-regulated DEGs were

Discussion
Leopard coral grouper is a valuable marine species in coral reef ecosystems, and an excellent research model for exploring its special skin color variations [5]. Fish skin coloration is a complex trait that is determined by genetic, cellular, physiological, and environmental factors [16]. The skin color in the fish, such as red leopard coral grouper and red tilapia, is the major factor of commercial value determination in Asian countries [45,46]. Therefore, it is essential to understand the genetic mechanism of color variations in these aquaculture species. To date, several investigations have reported many candidate genes associated with color changes through RNA-Seq in different non-model teleost fishes, such as common carp, red tilapia, and Taiwanese loach [41,42,47,48], as well as in model organisms such as zebrafish (Danio rerio) and medaka (Oryzias latipes) [49,50].   [42,45]. Previous studies reported erythrophore/xanthophore synthesis (slc7a11 and slc24a5), try, tryp1, dct, mitfa, and sox10 are well-recognized melanophoremarkers in red crucian carp and red tilapia [42,51].
Many SLC genes regulating the transport of substances on the cell membrane have been involved in skin pigmentation in humans and fish, such as slc45a2, slc24a5, and slc7a11. slc7a11 encodes the cystine/glutamate exchanger xCT that increases pheomelanin synthesis in red skin formation [42,51,52]. slc24a5 plays an essential role in regulating melanophore development in several vertebrates, such as zebrafish and common carp, so knockdown of the slc24a5 could block melanin synthesis to generate albino or golden phenotype [53,54]. slc45a2 participates in the membrane-associated transporter protein (MATP), which is likely involved in intracellular processing and trafficking of melanosomal proteins [55].  Table S1 Fig. 6 Comparison of gene expression patterns obtained using comparative transcriptome analysis and qRT-PCR. Nine genes were identified as DEGs between brown (control group) and red individuals (experimental group) (n = 4 for each group). Expression of target genes was normalized to β-actin as a reference gene. The Y-axis shows the relative mRNA expression levels (means ± SD). Statistically significant differences compared to control group are presented, with *: P < 0.05; **: P < 0.005; ***: P < 0.001 In the present study, we performed a comparative transcriptomic analysis between the red and brown skin of leopard coral grouper, and 1,089 significant DEGs were identified. Through comparative transcriptome analysis and qRT-PCR, we found that mRNA expression levels of epha2, sema6d, acsl4, slc7a5, hipk1, nol6, and timp2 were greater in the red individuals compared to brown individuals, and that mRNA expression levels of slc25a42 and kdf1 were lower in the red individuals compared to brown individuals. Two candidate SLC genes (slc7a5 and slc25a42) were identified as significant DEGs for skin color formation in leopard coral groupers. slc7a5 was an up-regulated DEG in red individuals, and slc25a42 was an up-regulated DEG in brown individuals. slc7a5 is an  Table 3 List of primers used for qRT-PCR. amino acid transporter that regulates the mTOR pathway, enhancing pheomelanin synthesis through activating xCT expression in mice [56]. The SLC25 family member is the largest solute transporter family in humans and transports solutes across the inner membrane of mitochondria [57]. slc25a38 contributed to red color development in spider mites [58], which may imply that slc25a42 was the potential pigmentation-related gene in leopard coral grouper. epha2 is a tyrosine kinase, which belongs to the family of Eph receptors. It was reported that epha2 inhibits MAPK and AKT pathways in human lens epithelial (HLE) cells [59], so the gene may prevent the MAPK pathway from activating the downstream gene of melanocyte inducing transcription factor (MITF), which is well known to participate in the enzymatic conversion of tyrosine to melanin [42,59]. acsl4 encodes a protein associated with lipid biosynthesis and fatty acid degradation [60]. Recent studies have indicated that acsl4 plays a critical role in activating ferroptosis that leads to cell death by iron-dependent lipid peroxidation, and also inhibiting melanin synthesis [61]. epha2 and acsl4 were up-regulated DEGs, which may inhibit melanin synthesis in red skin color individuals. In addition, melanin plays an essential role to protect animals from UV radiation and environmental challenges. Keratinocyte-derived factors are involved in regulating the proliferation and melanogenesis in mammals [62,63]. In our results, kdf1, a kind of keratinocyte-derived factor, had a higher expression in brown skin color individuals.
Taken together, most skin color-associated genes found in the study are involved in melanin synthesis which has been reported in previous studies. The genetic mechanisms of pigment cells are similar between fish and humans, and some skin color-associated genes in fish models, such as slc25a42 and MATP, have contributed to understanding the genetic mechanism in human skin [64].

Conclusions
This study conducted comparative transcriptome analysis in different skin colors of leopard coral groupers, and results showed that red individuals were different from brown individuals. The results of qRT-PCR showed that nine candidate genes were associated with the formation of skin coloration. To conclude, our results provided useful genetic resources for future studies regarding the genetic mechanisms of skin coloration in leopard coral grouper, and also to assist breeders to conduct molecular assisted selection in leopard coral grouper farming.

Sample collection
Fish were wild-caught from the exclusive economic zone of Penghu Sea in Taiwan. Skin samples were immediately immersed in the RNA keeper reagent (Protech, Taiwan), and further stored at -80℃. Six fish with different skin colors, including three red individuals and three brown individuals, were selected for RNA extraction. The average body weight of fish was 0.7 ± 0.6 kg, and average body length was 38.1 ± 9.3 cm.

RNA isolation
Total RNA was extracted from the fish skin tissue. The skin sample was homogenized in 1 mL of TRIzol reagent (Invitrogen, USA) containing stainless steel beads, following the manufacturer's protocol. The quality and quantity of total RNA were determined by a DS-11 spectrophotometer (DeNovix, USA), and run on 1% agarose gel electrophoresis. The RNA samples were stored at -80℃ for further experiments.

Library construction and RNA sequencing
RNA extracted from fish skin tissue was prepared for library construction. The RNA integrity was assessed by a Qsep400 (BiOptic, Taiwan) with quality number (RQN) values ≥ 7.0. A paired-end (PE) sequencing library was constructed using Illumina TruSeq RNA Library Prep Kit v2 (Illumina, USA), following the manufacturer's protocol. The obtained library was sequenced using the Illumina NovaSeq 6000 platform (Illumina, USA) with 2 × 150 bp PE reads.

Transcriptome assembly
To eliminate errors, adapters and low-quality reads from the raw data were removed by fastp [65]. The Phred score (Q20 and Q30) and GC content of the clean data were calculated. Sequence quality metrics were assessed using FastQC (http:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/) software. Subsequently, all clean reads were aligned to the reference P. leopardus genome [31] using STAR [66]. The coordinate sorted bam file was assembled using Trinity program [67] for reference-guided de novo assembly. While the assembled data contains abundant duplicate transcripts, CD-HIT-EST clustering application [68] was used to remove redundant transcripts with 95% identity, and establish the final assembled transcripts. Potential protein-coding genes were identified from assembled transcripts using Transdecoder (https:// github. com/ Trans Decod er/ Trans Decod er).

Differential expression analysis
Each PE read from the RNA-Seq was mapped with Bow-tie2 [70] to the final assembled transcripts, and the quantification of the read count was generated using RSEM [71]. The read counts were normalized by the TMM method implemented in the R package "edgeR" [72], and generated PCA and heatmap. Differential expression analysis was performed using "edgeR", and p-values were adjusted using the Benjamini-Hochberg (BH) method to control the FDR. Genes with FDR < 0.05 and an absolute value of log 2 FC ≥ 2 were considered to be DEGs, and were used in the subsequent analysis. The GO and KEGG enrichment analysis of DEGs were performed using KOBAS [73].

Quantitative real-time polymerase chain reaction (qRT-PCR)
cDNA synthesis was generated from total RNA using SuperScript ® IV First-Strand cDNA Synthesis Reaction (Invitrogen, USA), and random hexamer primers (Invitrogen, USA). Diluted cDNA (10 ng/µL) was used as a template for qRT-PCR. β-actin was used as the housekeeping gene for the normalization of gene expression levels. The primers used for target and reference genes were designed based on transcriptome sequence using Primer-BLAST [74] with an amplicon size of approximately 70-240 bp. qRT-PCR analysis was performed in the QuantStudio ™ 3 real-time PCR systems (Thermo Fisher Scientific, USA) and the reactions were carried out using PowerUp ™ SYBR ™ Green Master Mix (Thermo Fisher Scientific, USA). The conditions for reactions were 50℃ for 2 min, 95℃ for 2 min, followed by 40 cycles of 95℃ for 15 s, 60℃ for 15 min, 72℃ for 1 min. To ensure a single amplicon reaction, a melting curve analysis was performed at the end of each run. All the experiments were conducted with three biological replicates, and each group was conducted with four samples. The relative gene expression was calculated using the comparative threshold cycle (2 −△△CT ) method [75]. Student t-test was performed to determine the significant differences between red and brown groups. The Pearson correlation coefficient between RNA-Seq and qRT-PCR was calculated using "cor.test" function in R.